Pseudogap and Fermi-arc Evolution in the Phase-fluctuation Scenario 
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Pseudogap phenomena and the formation of Fermi arcs in underdoped cuprates are numerically 
studied in the presence of phase fluctuations that are simulated by an XY model. Most importantly 
the spectral function for each Monte Carlo sample is calculated directly and efficiently by the 
Chebyshev polynomials without having to diagonalize the fermion Hamiltonian, which enables us to 
handle a system large enough to achieve sufficient momentum/energy resolution. We find that the 
momentum dependence of the energy gap is identical to that of a pure d-wave superconductor well 
below the KT-transition temperature (Tkt ) , while displays an upturn deviation from cos k x — cos k y 
with increasing temperature. An abrupt onset of the Fermi arcs is observed above Tkt and the arc 
length exhibits a similar temperature dependence to the thermally activated vortex excitations. 
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One of the most important open issues in the research 
of high temperature superconducitity is the origin of the 
so-called pseudogap and its relation to the superconduct- 
ing gap. One viewpoint regards the pseudogap as pre- 
formed Cooper-pair state without long-range phase co- 
herence, yet another attributes its origin to a new phase 
of matter. The most direct evidence of the pseudogap 
in hole-doped cuprates comes from experimental tech- 
niques such as angle resolved photoemission spectroscopy 
(ARPES) and scanning tunneling spectroscopy (STS), 
which measure the low-energy single-particle spectrum 
and relate the loss of spectral weight to the pseudogap 
opening. [l|-|3j One peculiar property of the pseudogap 
phase revealed by ARPES is the truncated Fermi sur- 
faces termed as 'Fermi arcs'. [J] The Fermi arcs emerge 
abruptly around the nodal region above T c , and extend 
out to the Brillouin zone edge with increasing tempera- 
ture and form a complete Fermi surface for temperature 
above T* . [4- 6] This evolution of the Fermi surface shows 
a sharp distinction from that of conventional BCS super- 
conductors. 

The evolution of the Fermi surface in the pseudogap 
phase has been addressed [7-9] within the phase fluc- 
tuation context. The formation of Fermi arcs was at- 
tributed to the pile-up of low-energy spectral weight 
due to the Dopplcr effect of the fluctuating suppercur- 
rcnt. 7] However such analytic results are based on a 
semiclassical approximation [lfj which was argued to 
be uncontrolled 
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Numerical treatment of the ef- 
fect of phase fluctuations involved a combination of the 
Monte Carlo (MC) simulation of the Landau- Ginzburg 
free energy and the exact diagonalization of the fermion 
Hamiltonian. [8( However, the competition between the 
superconducting and antifcrromagnetic order parameters 
in the presence of quenched disorder was complicated and 
the momentum resolution was relative low due to the 
time-consuming numerical diagonalization procedure. 



In this work, we 1) examine the momentum and tem- 
perature dependence of the spectral gap and 2) tem- 
perature dependence of the Fermi-arc length within 
the phase-fluctuating scenario. The fluctuating phases 
are simulated phenomcnologically by a classical XY 
model. 12J LL2[ The XY free energy is 
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where the phase-stiffness constant J is the energy scale of 
the Kosterlitz-Thouless transition. tp\ denotes the phase 
of the superconducting order parameter defined on the 
lattice site i. S is the unit vector along x or y direction 
connecting i to its nearest neighbors (NN). The Hamil- 
tonian of non-interacting electrons moving in a spatially 
fluctuating pairing field as described by Eq. ([T]) reads 

6 = E(-^-^) c U-+E( Ai J c it c Ji+ iJ - c -)' 



U.o' 



>d 






y + M&j 



J )(2i 



where £y denotes hopping integral, and fj, represents the 
chemical potential. The pairing potential Ay is defined 
on the NN bond (i, j = i + 6). Here we fix the ampli- 
tude of the pairing potential while concentrate on the 
role played by the phase fluctuation, i.e A^i+a = Ae' 95 ' 1 ' 
with A kept constant. We use the expression 9x,y(i) = 
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to relate the phase of the pairing 



potential with i^j, which gives rise to pure d-wave pair- 
ing state in the low temperature limit when the thermal 
fluctuation is negligible small. 

In our numerical study, a set of phase configura- 
tions {fi} is first generated by the classical Monte 
Carlo simulation according to the probability distribu- 
tion exp(— F({(fi})/T). Then the single-electron spectral 



function of the BdG Hamiltonian ([2]) for each individ- 
ual phase configuration, which is denoted by A v {k., ui) in 
short, is calculated and at last averaged over the sam- 
pled phase space to obtain the expectation value of the 
spectral function A(k,uj) — (A v (k,w)). The bottleneck 
in this algorithm mainly lies in the numerical calculation 
of the spectrum of the BdG Hamiltonian ([2]). Previous 
studies [Si |l2| attempted to solve the problem by direct 
diagonalization of the BdG Hamiltonian, whose workload 
is 0(N 3 ) with N = L 2 the lattice size and rather time- 
consuming when L is large. In this paper, we deploy the 
Chebyshev expansion approach [14J to directly calculate 
the spectral function, which circumvents the difficulty in 
fully diagonalizing a large Hermitian matrix. As illus- 
trated in the following, the computational load for each 
A lf> (k, ui) is approximately O(MN) with M « N, which 
is a significant advance regarding to numerical studies of 
the above problem. 

The single-particle spectral function corresponding to 
the quadratic Hamiltonian of Eq. ([2J can be calculated 
from the imaginary part of the Green's function 
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A v {v,lo) = Im(i/| (w + ir\ — H v ) l \v). 
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Here more generally, | v) denotes any single-particle wave 
vector concerned with. \v) can be chosen as |k) = 
— ^^;e lkl |i) with k the momentum in the first Bril- 
louin zone, and accordingly the obtained ^4(k, u) are re- 
lated to the evolution of the Fermi surface as observed by 
ARPES. If \u) — |i), A(i, w) is actually the local density 
of state related to the STS measurements. According to 
the kernel polynomial method |14j], A v (v,u}) can be ex- 
panded in truncated series of the Chebyshev polynomials 
of the first kind T m , 
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where the expansion coefficients 
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H m = (v\T m (H / s)\v) 
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are Chebyshev moments. g m — sinh[A(l — 

m/M)]/ sinh(A) is the Lorentz kernel factor which im- 
proves the convergence of the truncated series and damps 
Gibbs oscillations. A good value of the kernel parameter 
is A = 4. The scaling factor s fits the spectrum of the 
Hamiltonian H/s into the interval [—1, 1]. We choose s 
to be a slightly larger than the bandwidth of the BdG 
Hamiltonian @. The energy resolution of the spectral 
function is estimated by e = As/A/. Most computational 
effort is spent in the calculation of the moments fi m ac- 
cording to Eq. ([5]), which reduces to matrix- vector multi- 
plications after taking advantage of the recursion relation 
T m {x) — 2xT m ^i(x) — T m _2(x). Considering that the 
BdG Hamiltonian ^ is sparse, the cost of matrix- vector 



multiplication is an order O(N) process and the calcula- 
tion of M moments requires only O(MN) computational 
operations, which even enables us to handle larger lat- 
tices using a desktop computer. Furthermore, recursive 
relations of Chebyshev polynomials T2,„ = 2X^ — 1 and 
T2m+i = 2T m T TO+ i — Xi enable us to obtain two moments 
per matrix- vector multiplication. 

We now present the results of our simulations. Typical 
calculation is done on a 96 x 96 lattice. Combined Monte 
Carlo method is used to generate equilibrium configura- 
tions of phase angles {ifii}. The KT-phase transition is 
observed with Tkt ~ J by studying the phase correla- 
tion and the vortex density as a function of temperature, 
which is consistent with our previous work. [9J For each 
temperature, the first 10 5 MC sweeps are dropped to 
equilibrate the system. 1000 phase configurations every 
1000 MC sweeps are used to get thermal averages, which 
are sufficient to reduce the statistical errors. For each 
phase configuration, the corresponding spectral function 
is calculated by Chebyshev expansions with the trunca- 
tion M = 4096, which results in an approximate energy 
resolution e/t ~ 0.004. Here t denotes the NN hop- 
ping integral. Most calculations are performed for the 
particle-hole symmetric case (t' = /j = 0) for simplicity. 
A = O.li which results in the superconducting gap 0.4i 
for the particle-hole symmetric case in the zero temper- 
ature limit. The resulting Fermi surface has a diamond 
shape, part of which is shown in the inset of Fig. [TJa). 
On the node-to-antinode segment of the Fermi surface, 
49 equally-spaced Fermi momenta are chosen to achieve 
sufficient momentum resolution. 

Figure [1] shows the energy distribution of the spec- 
tral function A(k, lu) for several Fermi momenta for il- 
lustration. Well below Tkt with T/J = 0.5 we find 
sharp spectral peaks, indicating well-defined Bogoliubov 
quasiparticle excitations. Furthermore the superconduct- 
ing gap vanishes only at the node, away from which 
the gap gradually opens and increases to its full size at 
the antinodal point. With increasing temperature from 
Fig. QJa)-(c), the spectral lines are broadened due to 
the strengthened scattering effect from the pair fluctu- 
ations. The energy gaps of k's near the antinode (see 
the evolution of the green and cyan curves) keep open- 
ing while that near the node (see red curves of (a)-(c) 
and (d)) closes as temperature rises across Tkt ~ J, 
indicating the formation of pseudogap at the antinodal 
region and the growth of Fermi arc near the nodal re- 
gion above Tkt- From the calculated spectral function 
in momentum- energy space, one can suggest two types 
of gap: the spectral gap A sg (k) which is half of the 
peak-peak separation [15j and the fitted gap A/ g (k) by 
numerical fitting with a simple phenomenological model 
A(k, to) = — Im[u) — £k — S(k, lo)]" 1 /tt with the modified 
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FIG. 1: Energy and momentum dependence of the spectral 
function (a), at T/J — 0.5 in the superconducting state, (b), 
at T/J — 1.0 around the KT phase transition and (c), at 
T/J — 1.6 in the pseudogap state. The inset of (a) shows 
the Fermi surface (black line) in first quarter of the first Bril- 
louin zone together with some selected Fermi momenta whose 
colors correspond to those of the curves in (a)-(c). (d), Varia- 
tion of the spectral function with temperature for fixed Fermi 
momentum (red dot in the inset) superimposed by their fits 
(black curves) according to Eq. ([6]) 



BCS self-energy, dHH | 
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where £k = — 2£(cos k x + cos k v ) — ^t 1 cos k x cos k y — n the 
normal-state electron dispersion relation. A/ 9 (k) and 
Tk are fitting parameters. Tk takes into account of the 
lifetime broadening due to the scattering effect from the 
pair fluctuations. Our model of a momentum-dependent 
scattering rate T^ is more general than that considered 
in Ref. [17L fl8j| . Fig.QJd) shows the typical fitting results 
together with the raw data calculated from Monte Carlo 
simulations. 

Figure [2] shows the dispersion of the spectral gap 
A ss (k) and fitted gap A/ S (k). At T/J = 0.5 (and 
below) both gap functions show simple d x 2_ y 2 form, 
| cos fez — cosky\/2, indicating that the effect of the pair 
fluctuation is rather weak well below Tkt- With in- 
creasing temperature until T/J =1.0 which is approxi- 
mately equal to Tkt, both gap functions deviate from the 
pure d-wave form and exhibit an upturn behavior due to 
increased spin-wave-type superconducting fluctuations. 
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FIG. 2: Momentum dependence of energy gap at four dif- 
ferent temperatures T/J — 0.5, 1.0, 1.1, 1.2(black dot, red 
square, blue diamond, and green triangle), (a) spectral gap 
(half the peak-to-peak distance) plotted as a function of <f), 
i.e. the angle of k relative to the nodal point, (b) fitted gap 
(see Eq.[SJ) plotted as a function of | cos(fc;c) — coa(k y )\/2. The 
gap is normalized to its value at the antinode. Inset of (b) 
shows the angle dependence of T k (see Eq. [5| 



This deviation increases more rapidly with temperature 
from T/J = 1.0 to 1.2, which indicates that the vortex- 
like excitations are more efficient in Doppler-shifting the 
single-electron spectral functions. From Fig.[2jb), the be- 
havior of Af g changes gradually with temperature. How- 
ever by examining A sg as in (a), we find that an emergent 
"arc", i.e. a Fermi segment with zero spectral gap, arises 
abruptly around Tkt and grows rapidly with tempera- 
ture from T/J — 1.0 to 1.2. This expansion of a Fermi 
node (below Tkt) to an Fermi arc (above Tkt) can be 
attributed to thermal proliferation of vortices. 

The length of Fermi arc can be quantified by differ- 
ent ways. Experimentally one can locate the specific k 
where the spectral gap begins to open by examining the 
k dependence of A sg (k) , or by examining the loss of spec- 
tral weight. [6j Here we focus on the fact that Fermi arcs 
are segments of Fermi surface where the single-electron 
spectral function at the Fermi energy has non-vanishing 
value. Therefore A(k, u = 0) is studied with k along 
the underlying Fermi surface of the non-interacting elec- 
trons. The results are shown in Fig. [3] At T/J = 0.5 
we find a sharp peak around the nodal point, which de- 
cays rapidly along the node-to-antinode direction indi- 
cating point-like Fermi surface in the superconducting 
state. This peak is broadened with increasing tempera- 
ture and accordingly phase fluctuations. Compared with 
the Gaussian- type decay obtained in Ref. [7| , we find that 
A(k, 0) can be better fitted by a Lorentzian-type func- 
tion / 2 as shown in Fig. GHa), where (f>u is the angle 
of k relative to the node, and a, <f) c are fitting parame- 
ters. This is not accidental but compatible with Eq. ([6]). 
For k along the Fermi surface and near the gap node, 
A^ ff (k) m (\/2tt — &_F)uA0k an d after considering that 
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FIG. 3: (a), Momentum dependence of the spectral function 
at the Fermi level A(k.,w — 0) with k along the underlying 
Fermi surface at four temperatures T/J — 0.5, 1.0, 1.2, 1.6 
(black circle, red square, blue diamond, green triangle). <f>k 
denotes the angle of k relative to the node point (n/2, n/2). 
Solid colored lines are fits to a Lorentzian-type function ,..° j9 

(b), the characteristic width <j> c , as a function of temperature 
from T/J = 0.1 to 3.0. The inset shows the temperature 
dependence of vortex density. 
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Here 4> c can be regarded as the characteristic width of the 
decay of the spectral function and accordingly as a mea- 
sure of the length of Fermi arc. Fig. [SJb) shows the tem- 
perature dependence of (f> c . Below Tkt the arc length is 
negligible small, showing that the Fermi node is protected 
against weak phase fluctuations. Near T/J =1.0 there 
is an apparent jump of <p Cl indicating that the Fermi arc 
emerges immediately above Tkt together with the loss 
of phase coherence, which is consistent with experimental 
observation. [5] The relation of Fermi arc formation to the 
vortex proliferation is illustrated by the inset of Fig. EJb) 
where the temperature dependence of the vortex density 
is shown. The formation of vortices increases rapidly for 
temperatures near and above Tkt and exhibits a similar 
temperature behavior to the arc length indicating that 
vortex-type phase fluctuations are responsible for the for- 
mation and evolution of the Fermi surface. This result is 
ualitatively consistent with an earlier theoretical study 
in which the arc length is related to the characteris- 
tic width of the supercurrent distribution. However, our 
work has taken both the Doppler effect of the whirling 
supercurrent and the scattering effect of vortex centers 
into full consideration. 

In Fig. [4] we plot the intensity of the spectral function 
at zero energy as a function of k x and k y in the first quar- 
ter of the first Brillouin zone for a realistic hopping pa- 
rameter with t' /t = —0.3. Well below the KT-transition 
temperature with ignorable phase fluctuations, the Fermi 
surface is point like as in ideal d-wave superconductors. 




FIG. 4: The temperature dependence of the spectral function 
at the Fermi level, A(k,uj — 0) normalized to its value at the 
node point for realistic band parameters t' /t = —0.3. Here we 
only show the results in the first quarter of the first Brillouin 
zone. From (a) to (d), T/J = 0.5, 1.0, 1.6, 2.0. 



With increasing temperature we find arc-like Fermi sur- 
face with monotonically increasing arc length, which is 
qualitatively consistent with ARPES observation [g, [6j . 

In summary, the influence of thermal phase fluctuation 
on the single-particle spectral function has been studied 
via numerical techniques which combine the Monte Carlo 
method with the Chebyshev expansion approach. The 
dispersion of the spectral gap deviates from the simple 
d-wave form especially for T > Tkt when vortex-like 
fluctuations are dominant. Apparent onset of Fermi arcs 
are found near and above Tkt, with the arc length rises 
monotonically with temperature and displays close rele- 
vance to the vortex density. Our results give qualitative 
understanding of the ARPES observations of the pseu- 
dogap phase and Fermi arcs in cuprates. Quantitatively, 
specific temperature dependence of the phase stiffness J 
is also crucial and its determination may calls for start- 
ing from a microscopic Hamiltonian. This issue will be 
addressed in the future work. 
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